Resistive heater hysteresis control¶
In this section, the hysteresis control for resistive heaters is modelled. The resistive heater comprises an on/off controlled heating system, which means the resistive heater is either switched off or operates at its rated power. The COP of resistive heaters is 1, which means the amount of heat gain from the resistive heater is equivalent to its rated power $P_{heater}$ (e.g., 2kW).
The power of heat loss $\frac{d Q_{loss}}{dt}$ (unit: W) can be calculated as follows[1]:
\begin{align*} \frac{d Q_{loss}}{dt} = K \cdot (T_{in} - T_{out}) \end{align*}
where $\frac{d Q_{loss}}{dt}$ (unit: W) is the power of heat loss of the building. $T_{in}$ and $T_{out}$ represents the indoor and outdoor temperatures. $K$ (unit: W/K) denotes the building loss coefficient, which is the product between the thermal transmittance U-value (unit: W/($m^{2} \cdot K$)) and the surface area A (unit: $m^{2}$). The available U-value for each building component can be found in [2][3].
Finally, the indoor temperature changes can be calculated as follows (no internal gains) [4][5][6][7][8]:
\begin{align*} \frac{d T_{in}}{dt} = \frac{1}{m_{air}^{in} \cdot c_{air}} (P_{heat} - \frac{d Q_{loss}}{dt}) \end{align*}
where $\frac{d T_{in}}{dt}$ (unit: \degree C/s) shows the indoor temperature change. The building's air mass is denoted as $m_{air}^{in}$ (unit: kg).
References¶
[1]: Rogeau, A., Vieubled, J., Ruche, L. and Girard, R., 2024. A generic methodology for mapping the performance of various heat pumps configurations considering part-load behavior. Energy and Buildings, 318, p.114490.
[2]: Clauß, J., Stinner, S., Sartori, I. and Georges, L., 2019. Predictive rule-based control to activate the energy flexibility of Norwegian residential buildings: Case of an air-source heat pump and direct electric heating. Applied Energy, 237, pp.500-518.
[3]: Alimohammadisagvand, B., Jokisalo, J. and Sirén, K., 2018. Comparison of four rule-based demand response control algorithms in an electrically and heat pump-heated residential building. Applied Energy, 209, pp.167-179.
[4]: Saleh, A. and Mosa, M., 2016. Analysis of control strategies and simulation of heating systems using Simulink/Matlab potential. Journal of Thermal Engineering, 2(5), pp.921-927.
[5]: Xiang, S., Chang, L., Cao, B., He, Y. and Zhang, C., 2019. A novel domestic electric water heater control method. IEEE Transactions on Smart Grid, 11(4), pp.3246-3256.
[6]: Buonomano, A. and Palombo, A., 2014. Building energy performance analysis by an in-house developed dynamic simulation code: An investigation for different case studies. Applied Energy, 113, pp.788-807.
[7]: Tindale, A., 1993. Third-order lumped-parameter simulation method. Building Services Engineering Research and Technology, 14(3), pp.87-97.
[8]: Yao, R. and Steemers, K., 2005. A method of formulating energy load profile for domestic buildings in the UK. Energy and buildings, 37(6), pp.663-671.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.offline as pyo
import plotly.graph_objects as go
import matplotlib.dates as mdates
from plotly.subplots import make_subplots
import plotly.io as pio
pio.renderers.default = 'notebook_connected'
from datetime import timedelta
df = pd.read_csv('simulation_results_finished.csv')
df['Date'] = pd.to_datetime(df['Date']) # Convert 'Date' column to datetime
df.set_index('Date', inplace=True)
outdoor_temp_series = df['Dry-bulb temperature (°C)']
outdoor_temp_series = outdoor_temp_series.to_frame().reset_index()
# Convert time to seconds for simulation alignment instead of half-hourly intervals
outdoor_temp_series["time_seconds"] = (outdoor_temp_series["Date"] - outdoor_temp_series["Date"].iloc[0]).dt.total_seconds()
# Define constants
heater_power = 2000 # Heater power in watts (W)
room_volume = 180*2.6 # Room volume in cubic meters (m³)
air_density = 1.225 # Air density in kg/m³
specific_heat_air = 1005 # Specific heat of air in J/(kg·K)
heat_loss_rate = 30 # Heat loss rate to environment in W/K (depends on insulation)
initial_indoor_temp = 20 # Initial indoor temperature in °C
thermostat_temp = 20 # Thermostat set temperature in °C
hysteresis = 2.0 # Thermostat hysteresis in °C
simulation_time = 365 * 24 * 3600 # Simulation time for one year (in seconds)
time_step = 1 # Time step in seconds
air_mass_flow_rate = 0.1 # Air mass flow rate in kg/s
# Derived constants
air_mass = room_volume * air_density # Total air mass in the room (kg)
# Simulation variables
time = np.arange(0, simulation_time, time_step) # Time array for one year
indoor_temp = np.zeros_like(time, dtype=float) # indoor temperature array
heater_status = np.zeros_like(time, dtype=int) # Heater status array (1=ON, 0=OFF)
heater_temp = np.zeros_like(time, dtype=float) # Heater temperature array
heat_loss = np.zeros_like(time, dtype=float) # Heat loss array
# Interpolate outdoor temperature for simulation time steps
outdoor_temp = np.interp(time, outdoor_temp_series["time_seconds"], outdoor_temp_series['Dry-bulb temperature (°C)'])
# Initialize indoor temperature
indoor_temp[0] = initial_indoor_temp
# Simulation loop
for t in range(1, len(time)):
# Calculate heat losses
heat_loss[t-1] = heat_loss_rate * (indoor_temp[t-1] - outdoor_temp[t-1])
# Determine heater status based on thermostat control
if indoor_temp[t-1] < thermostat_temp - hysteresis / 2:
heater_status[t] = 1 # Heater ON
elif indoor_temp[t-1] > thermostat_temp + hysteresis / 2:
heater_status[t] = 0 # Heater OFF
else:
heater_status[t] = heater_status[t-1] # Maintain previous state
# Calculate heat gain from heater
heat_gain = heater_status[t] * heater_power
# Calculate heater temperature
heater_temp[t-1] = indoor_temp[t-1] + heat_gain / (air_mass_flow_rate * specific_heat_air)
# Calculate net heat change
net_heat = heat_gain - heat_loss[t-1]
# Update indoor temperature
temp_change = (net_heat * time_step) / (air_mass * specific_heat_air)
indoor_temp[t] = indoor_temp[t-1] + temp_change
# Save the results in a DataFrame
sim_results = pd.DataFrame({
'Time': time,
'Outdoor Temperature (°C)': outdoor_temp,
'Indoor Temperature (°C)': indoor_temp,
'Heater Temperature (°C)': heater_temp,
'Heat Gain Power (W)': heater_status * heater_power,
'Heat Gain Energy (kWh)': (heater_status * heater_power * time_step) / 3.6e6,
'Heat Pump Status': heater_status,
'Heat Loss Power (W)': heat_loss,
'Heat Loss Energy (kWh)': (heat_loss * time_step) / 3.6e6,
})
sim_results.set_index('Time', inplace=True)
sim_results.index = pd.to_datetime(sim_results.index, unit='s', origin='unix').strftime('%m-%d %H:%M:%S')
sim_results.to_csv('../Rule based control - heat pumps/dataset/resistive_heater.csv', index=True)
Visualisation¶
df = pd.read_csv('../Rule based control - heat pumps/dataset/resistive_heater.csv')
df.set_index('Time', inplace=True)
df.head()
| Outdoor Temperature (°C) | Indoor Temperature (°C) | Heater Temperature (°C) | Heat Gain Power (W) | Heat Gain Energy (kWh) | Heat Pump Status | Heat Loss Power (W) | Heat Loss Energy (kWh) | |
|---|---|---|---|---|---|---|---|---|
| Time | ||||||||
| 01-01 00:00:00 | 5.5 | 20.000000 | 20.000000 | 0 | 0.0 | 0 | 435.000000 | 0.000121 |
| 01-01 00:00:01 | 5.5 | 19.999245 | 19.999245 | 0 | 0.0 | 0 | 434.977350 | 0.000121 |
| 01-01 00:00:02 | 5.5 | 19.998490 | 19.998490 | 0 | 0.0 | 0 | 434.954702 | 0.000121 |
| 01-01 00:00:03 | 5.5 | 19.997735 | 19.997735 | 0 | 0.0 | 0 | 434.932054 | 0.000121 |
| 01-01 00:00:04 | 5.5 | 19.996980 | 19.996980 | 0 | 0.0 | 0 | 434.909408 | 0.000121 |
# Calculate heater runtime and energy usage
heater_runtime = np.sum(df['Heat Pump Status']) * time_step / 3600 # Hours
energy_usage = np.sum(df['Heat Gain Energy (kWh)']) # kWh
# Display results
print(f"Heater runtime: {heater_runtime:.2f} hours")
print(f"Total energy usage: {energy_usage:.2f} kWh")
Heater runtime: 1562.87 hours Total energy usage: 3125.74 kWh
# Function to filter data for a specific day
def filter_data_for_day(df, day):
start_time = pd.to_datetime(day)
end_time = start_time + pd.Timedelta(days=1)
return df[(df.index >= start_time.strftime('%m-%d %H:%M:%S')) & (df.index < end_time.strftime('%m-%d %H:%M:%S'))]
# Choose a specific day to display
specific_day = '2023-08-20 00:00:00' # Change this to the desired date format 'YYYY-MM-DD HH:MM:SS'
daily_profile = filter_data_for_day(df, specific_day)
# Create figure with secondary y-axis
fig = make_subplots(specs=[[{"secondary_y": True}]])
# Add traces
fig.add_trace(
go.Scatter(x=daily_profile.index, y=daily_profile['Indoor Temperature (°C)'], name='Indoor Temperature (°C)', line=dict(color='blue')),
secondary_y=False,
)
fig.add_trace(
go.Scatter(x=daily_profile.index, y=daily_profile['Outdoor Temperature (°C)'], name='Outdoor Temperature (°C)', line=dict(color='orange')),
secondary_y=False,
)
fig.add_trace(
go.Scatter(x=daily_profile.index, y=daily_profile['Heat Gain Power (W)'], name='Heat Gain Power (W)', line=dict(color='green', dash='dot')),
secondary_y=True,
)
# Add thermostat bounds
fig.add_trace(go.Scatter(x=daily_profile.index, y=[thermostat_temp]*len(daily_profile.index), line=dict(color='red', dash='dash'), name="Thermostat Temperature (°C)"), secondary_y=False)
fig.add_trace(go.Scatter(x=daily_profile.index, y=[thermostat_temp - hysteresis / 2]*len(daily_profile.index), line=dict(color='grey', dash='dash'), name="Thermostat Lower Limit (°C)"), secondary_y=False)
fig.add_trace(go.Scatter(x=daily_profile.index, y=[thermostat_temp + hysteresis / 2]*len(daily_profile.index), line=dict(color='grey', dash='dash'), name="Thermostat Upper Limit (°C)"), secondary_y=False)
# Add figure title and axis labels
fig.update_layout(
title_text=f"Indoor temperature, outdoor temperature, and heat gain power from {pd.to_datetime(specific_day).strftime('%m-%d %H:%M:%S')} to {(pd.to_datetime(specific_day) + timedelta(days=1)).strftime('%m-%d %H:%M:%S')}",
width=1000, # Set the width of the figure
height=700, # Set the height of the figure
legend=dict(
orientation="h", # Horizontal orientation
x=0.5, # x position of the legend
y=-0.4, # y position of the legend
xanchor="center", # Anchor the legend at the center
yanchor="bottom", # Anchor the legend at the bottom
bgcolor='rgba(255, 255, 255, 0.5)', # Background color of the legend
bordercolor='Black', # Border color of the legend
borderwidth=1, # Border width of the legend
itemwidth=30, # Width of each legend item
traceorder="normal" # Order of legend items
),
xaxis=dict(
tickformat='%H:%M', # Format ticks as hours and minutes
dtick=3600, # Set tick step to 1 hour (3600000 milliseconds)
tickangle=45 # Rotate x-axis labels for better readability
),
)
# Set x-axis title
fig.update_xaxes(title_text="Time (s)")
# Update y-axes
fig.update_yaxes(title_text='Temperature (°C)', secondary_y=False)
fig.update_yaxes(title_text='Heat Gain Power (W)', secondary_y=True, tickfont=dict(color='green'), titlefont=dict(color='green'))
# Show plot
fig.show()